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FOREWORD 


This  report  was  prepared  by  the  University  of  Dayton  Research  Institute  under  Air  Force 
Contract  No.  F33615-95-D-5029,  Delivery  Order  No.  0006.  The  work  was  administered  under 
the  direction  of  the  Nonmetallic  Materials  Division,  Materials  and  Manufacturing  Directorate, 
Air  Force  Research  Laboratory,  Air  Force  Materiel  Command,  with  Dr.  L.  Scott  Theibert 
(AFRL/MLBC)  as  Project  Engineer. 

This  report  was  submitted  in  October  2000  and  covers  work  conducted  from 
15  September  1999  through  31  July  2000. 


EXECUTIVE  SUMMARY 


A  three-dimensional  model  for  a  stress  analysis  of  woven  fabric  composites,  which  was 
derived  previously  based  on  Reissner’s  mixed  variational  principle,  was  solved  numerically  with 
a  finite  element  approach.  Since  the  mixed  model  calculates  the  stress  field  by  taking  variations 
of  displacement  and  stresses  independently  and  satisfying  equilibrium  of  stresses  pointwise, 
accurate  interlaminar  stresses  are  predicted  at  the  yarn  interface.  The  interfacial  continuity 
conditions  are  implemented  through  a  penalty  method  by  adding  an  additional  variational  energy 
of  two  constraint  conditions:  the  displacements  must  be  continuous  along  the  interface  between 
two  stacked  subregions,  and  interfacial  normal  and  shear  stresses  must  be  in  equilibrium  at  the 
interface. 

After  performing  the  thickness  integration,  the  three-dimensional  variational  energy 
equation  is  evaluated  for  each  yarn  (subregion)  two-dimensionally  with  16  stress-related  and  13 
displacement-related  unknown  variables.  Using  the  Rayleigh-Ritz  approximation  yields  a 
system  of  linear  equations  by  taking  derivatives  of  the  variational  energy  equation  with  respect  to 
the  independent  unknown  variables.  The  present  mixed  method  is  applied  to  the  analysis  of  both 
flat  and  woven  laminated  composites.  The  displacement  and  stress  results  of  the  present  method 
are  compared  and  validated  with  the  conventional  displacement-based  finite  element  solutions 
and/or  the  existing  analytic  solution. 

The  present  model  is  also  applied  to  the  analysis  of  stiffness  and  strength  of  flat  and 
woven  fabric  composites.  The  model  calculates  three-dimensional  effective  elastic  moduli  and 
predicts  failure  strengths  and  damage  modes.  The  failure  analysis  includes  residual  stress 
calculation  to  consider  the  hygrothermal  effect.  The  numerical  calculations  show  good 
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agreement  with  existing  experimental  and  numerical  results  on  both  flat  and  woven  laminated 
composites  with  various  yarn-waviness  ratios. 

The  emerging  ultra-lightweight  material,  carbon  foam,  was  modeled  with  the  three- 
dimensional  microstructures  to  develop  a  basic  understanding  of  the  performance  of  open-cell 
foam  materials.  The  model  can  describe  the  deformation  behavior  accurately  and  will  be  used  to 
investigate  the  failure  mechanism  of  the  cell  ligaments. 

Because  of  the  randomness  and  complexity  of  the  microstructure  of  the  carbon  foam,  the 
representative  cell  ligaments  are  first  characterized  in  detail  at  the  micro  structural  level.  The 
microstructural  characterization  will  then  be  correlated  with  the  macroscopic  bulk  properties  by  a 
statistical  approach.  A  series  of  databases  will  be  collected  for  various  size  and  spatial 
orientation  of  the  cell  ligaments,  as  well  as  the  property  variation  due  to  the  graphic  alignment 
along  the  longitudinal  direction  of  the  ligaments. 

Because  of  slenderness,  each  ligament  can  be  considered  as  a  beam,  and  the  tetrahedral 
cell  microstructure  with  four  ligaments  as  a  frame  structure.  The  four  beams  are  located  in  three- 
dimensional  space  under  arbitrary  loading  conditions.  The  cross  section  of  the  beam  varies  in 
size  and  material  properties  in  the  longitudinal  direction  along  the  ligament. 

Based  on  the  literature  review,  the  research  approach  of  modeling  the  carbon  foam 
blowing  process  was  developed.  The  model  consists  of  three  concepts:  (1)  nucleation  of 
microcellular  that  determines  the  relationship  of  the  cell  number,  gas  kind,  temperature  and 
pressure,  (2)  bubble  growth  that  calculates  bubble  dynamic  size  and  shape  so  the  relationship  of 
the  process  parameters  and  foam  properties  can  be  determined,  and  (3)  mold  filling  that 
simulates  the  process  of  the  foam  filling  a  mold  cavity.  A  finite  element  code  is  being  developed. 
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1.  THREE-DIMENSIONAL  MODELING  OF  WOVEN  COMPOSITES 


1.1  INTRODUCTION 

To  achieve  the  optimum  structural  properties  of  state-of-the-art  fabric  reinforcements  of 
woven  composites,  there  is  a  need  to  develop  a  basic  understanding  of  deformation  and  damage 
mechanisms.  A  model  for  three-dimensional  stress  analysis  of  woven  fabric  composites  has 
been  proposed  by  Roy  [1]  to  obtain  reliable  three-dimensional  stress  fields,  especially 
interlaminar  stresses  along  interfaces  between  yarns.  The  model  is  formulated  based  on 
Reissner’s  mixed  variational  principle  to  take  independent  variations  on  the  stress  and  the 
displacement  components  [2-4].  The  in-plane  stresses  within  a  yarn  are  assumed  to  vary  linearly 
in  the  thickness  direction,  and  the  expressions  for  the  interlaminar  stresses  are  obtained  by 
satisfying  the  three-dimensional  equilibrium  equations.  After  performing  the  thickness 
integration,  the  three-dimensional  variational  energy  equation  becomes  a  two-dimensional 
equation.  The  variational  equation  is  expressed  with  16  stress-related  and  13  displacement- 
related  unknown  variables.  In  this  model,  an  accurate  calculation  of  the  interlaminar  stresses  at 
the  yarn  (subregion)  interface  can  be  achieved  (except  near  the  point  of  singularity)  by  satisfying 
the  interfacial  traction  continuity  conditions  and  the  equilibrium  of  stresses  pointwise. 

Our  present  work  establishes  a  mixed  finite  element  analysis  (FEA)  based  on  the  mixed 
variational  principle.  Total  variational  energy  is  obtained  by  accumulating  the  energy  for  all 
yam  and  matrix  subregions.  The  subregions  are  further  discretized  into  finite  elements  in  a  plane 
perpendicular  to  the  thickness  direction.  The  interfacial  continuity  conditions  are  implemented 
through  a  penalty  method  by  adding  an  additional  variational  energy  of  two  constraint 
conditions:  the  displacements  must  be  continuous  along  the  interface  between  two  stacked 
subregions,  and  interfacial  normal  and  shear  stresses  must  be  in  equilibrium  at  the  interface. 
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Two  large  numbers  of  penalty  parameters  enforcing  the  displacement  and  stress  continuity  are 
employed  carefully  to  avoid  numerical  errors.  The  solution  of  the  variational  energy  equation  is 
obtained  by  using  the  Rayleigh-Ritz  approximation  with  polynomial  shape  functions. 

The  present  mixed  method  is  applied  to  analyze  a  flat  laminated  composite  with  a  free 
edge  and  a  representative  volume  element  (RVE)  of  plain-weave  composites.  The  displacement 
and  stress  results  of  the  present  method  are  compared  with  the  conventional  displacement-based 
finite  element  solutions  and/or  the  existing  analytic  solutions.  The  reliable  prediction  of  the 
stress  field  by  the  present  method  is  used  to  calculate  stiffness  and  strength  of  the  flat/woven 
laminated  composites.  Three-dimensional  effective  elastic  moduli  are  calculated  for  several 
flat/woven  laminated  composites  and  compared  with  existing  experimental/numerical  results. 
Meanwhile,  a  discrete  damage  analysis  is  achieved  to  calculate  first-ply  and  last-ply  failure  loads 
as  well  as  their  damage  modes.  Failure  strengths  are  predicted  by  considering  not  only 
mechanical  stresses  but  also  residual  stresses  that  are  significantly  influenced  by  hygrothermal 
effects.  The  numerical  predictions  on  the  failure  strength  and  the  damage  mode  are  compared 
with  experimental  results  that  were  previously  observed  on  flat  laminated  composites  and  woven 
model  laminates  with  one-dimensional  yarn  crimping. 

1.2  FORMULATION  OF  THREE-DIMENSIONAL  MODEL 

1.2.1  Modified  Variational  Energy  Equation 

The  variational  energy  equation  evaluated  for  a  given  (fc-th)  subregion  is  written 
as 


+  F3v  +  F4v' 


+  F5w  +  F6w 


+  F7  vv)]  1  dxdy 


*y 

+  ||  [(^52  ~KxPl2  -h2.yP62)U2-{p5l  ~hl,xPll  ~hl,yP6l)Ul 

xy 


(1) 
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+  (p42  ~KxP62  ~h2,yP 22)^2  ~(p41  ~K,XP61  -kl,yP2l)vi 
+  (P32  ~  h2,xP52  ~  h2,yP42  )  W2  ~  (p3l  ~  KxP  51  ~  KyP4l  )  Wl  dxdV 

+  j  \P6F  +  P62U  ^  +  P2,V  +  P22V ‘  +  P4,W  +  P 42  ^  +  P43™)(h2  ~  kl  Hi  2  dX 

x 

+  |  {(pijU  +  Pnu*  +  PeiV+  P62V  +  pnw  +  p52w*  +  p53w)(h2  -hj  )]a’|  dy 

y  x‘ 

-  j  [fe2M2  +%2*2  +  f,2W2  )-  fe/M/  +^v/V/  +  ^7  f"  ^  =  0 

A 

where  jll;  and  x,;  are  defined  in  Roy  [1]. 

The  interfacial  continuity  condition  dictates  that  the  displacements  must  be 
continuous  along  the  interface  between  two  stacked  subregions  (k-th  and  Z-th  subregions),  and 
interfacial  normal  and  shear  stresses  must  be  in  equilibrium,  as  Figure  1  shows. 


Figure  1.  Displacement  and  Stress  Continuity  at  the  Interface  between  k- th  and  /- th  Subregions. 


By  setting  the  interfacial  normal  stress  as  <$  3  and  interfacial  shear  stresses  as  6  4  and  c?  5 ,  the 
interfacial  continuity  condition  provides  the  following  constraint  conditions: 

(1)  Displacement  continuity: 

u?-u?=0 

V™-V?=0  (2) 
W*'-w?=0 
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(2)  Normal  and  shear  stress  continuity: 


df  -df  =0 


drt)  -d 

'“’4 


-^5 


(/)  =0 
(/)  _0 


(3) 


where  <y:  and  <x  '  are  the  interfacial  stress  components  at  the  k-th  and  /-th  subregions, 
respectively.  Note  that  the  interfacial  stresses  are  evaluated  in  a  local  coordinate  system  whose 
planar  coordinates  are  parallel  to  the  interfacial  surfaces.  These  local  stress  components  along 
the  interfacial  surfaces  are  related  with  stress  components  in  the  global  coordinates  by  the  slopes 


of  the  interfacial  surfaces  in  x-  and  y-directions  ( h\  ’x  and  ti2  ’ ).  Stress  transformation  using  the 
direction  cosines  of  the  interfacial  surface  vectors  yields  the  following  stress  constraint  equations 
in  the  global  coordinate  system: 

d<*>-of=0  =>  -<s?)  =  0 


a?'-$n=0  =>  a 


(*)  _(/)  lik)  /  (&)  _(/)  \  jAk) 

4  ^4  ^2  x  *  ^6  )  ^2,  v  *  \^2  ^2  )  ^ 


(4) 


df-df=0  =*  =  o 


To  impose  the  constraint  conditions  for  displacement  and  stress  continuity,  one 
can  substitute  them  into  equation  (1)  directly  to  formulate  an  irreducible  form.  It  is 
straightforward  to  substitute  the  displacement  continuity.  However,  it  turns  out  that  substituting 
the  stress  continuity  requires  an  extremely  involved  algebraic  manipulation.  Moreover,  when 
obtaining  numerical  solutions  by  using  polynomial  shape  functions,  the  restriction  of  excessive 
continuity  for  stresses  should  be  avoided  at  singularities  and/or  at  abrupt  changes  in  material 
properties.  The  imposition  of  such  continuity  would  likely  produce  erroneous  and  usually  highly 
oscillating  results  [5]. 
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Instead  of  using  the  irreducible  form,  a  penalty  approach  is  employed  by  adding  a 
new  energy  term  (  Jc )  for  the  constraint  conditions,  with  penalty  parameters  ( OC,  and  oc, ),  which 
yields  the  following  modified  variational  energy  equation: 

J  =  tr“  +  L7"’  (5) 

k  =1  k=l 

where 

/<*>=-• ay  ■{{  [(u(k)  -«y>)2  +  (vf  -v\l))2+{w\k)-w^)2]dxdy 

(6) 

+  ya y\\\6f-6ff  +  {6r-6f)’  +  (6r-6ff\dxdy 

xy 

and  M  is  the  number  of  subregions  in  the  thickness  (z.)  direction.  Two  large  numbers  of  (X,  and 
OC2  enforce  the  displacement  and  stress  continuity,  respectively.  However,  a2  must  be  selected 
carefully  to  avoid  the  excessive  continuity  for  stresses.  Because  of  the  nature  of  the  mixed 
formulation  for  the  stress  and  the  displacements,  erroneous  results  in  stress  may  ruin  the  ones  in 
displacement,  and  vice  versa.  The  effect  of  the  penalty  parameters  will  be  discussed  later. 

Because  of  the  complexity  of  the  modified  variational  equation,  it  is  more 
desirable  to  obtain  the  solution  numerically  rather  than  analytically.  Using  the  Rayleigh-Ritz 
approximation  can  yield  a  system  of  linear  equations  that  is  solvable  numerically.  There  are  two 
possible  approaches,  finite  element  or  finite  difference,  which  can  be  taken  to  solve  the  system  of 
equations  numerically,  and  the  former  is  taken  in  this  study. 

Because  of  the  through-the-thickness  (z)  integration  during  formulation,  the 
modified  mixed  variational  equation,  equation  (5),  is  only  a  function  of  x  and  y,  and  so  are  the 
29*Ns  unknown  variables  ( C.k)(x,  y) ,  i  =  1 , •  •  •,  29 )  for  the  k-th  subregion,  where  Ns  is  the 


7 


number  of  subregions.  Among  29  unknown  variables  for  each  subregion,  16  are  for  the  stress 
components,  and  13  are  for  the  displacement  components,  as  in  equation  (8).  The  variational 
equation  is  then  discretized  in  x-  and  y-directions  for  the  finite  element  formulation.  The 
unknown  variables  are  collected  in  a  vector,  as  follows: 

{c,l»(x,y)}=|P”|  (7) 

where 

P  ~  {P 11  ’  P 12  ’  Pll  ’  P22  ’  P  3  V  P  32  ’  P  33’  P  34’  P4I  ’  P 42  ’  P  43  ’  P 51  ’  P  52’  P53’  P  61’  P  62  } 

and  (8) 

j(£)  f —  »  —  *  —  *  ~  \k)T 

a  —\u,u  ,  Uj,  u2,  V,  V  ,  Vj,V2,  W,W  ,W,Wj,  w2f 

Each  of  the  unknown  variables,  C\k)(x,y) ,  are  then  interpolated  with  their  nodal 
contribution,  ,  by  shape  functions,  as  follows: 

Nen 

Pik) (x,  y)  =  Yj  P IP  ■  NPj (x, y) 
j=i 

(9) 

Nen 

dy(x,y)  =  '^d,i'-Ndj(x,y) 
j=l 

where  Nen  is  the  number  of  nodal  points  in  an  element,  and  N  pj  and  Ndj  are  the  shape  functions 
for  the  stress  and  displacement  degrees  of  freedom,  respectively.  The  shape  functions  can  be 
chosen  as  the  linear  polynomial  for  4-node  quadrilateral  elements  ( N en  =4),  quadratic 
polynomial  for  8-node  serendipity  elements  ( Nen  = 8 ),  etc. 

The  nodal  values  of  the  unknown  variables  for  each  finite  element  are  collected  in 
a  vector,  as  follows: 
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(10) 


C(i)(x,y)  =  {c^>}=4)1,d1,p2,d2,-,pj,dj,-,pNen,dNen}w 

where 

P  j  —  {p  11  ’  P 12  ’  ^27  ’  P 22’  P 31  ’  P 32  ’  ’  P 34  ’  P 41  ’  P 43’  P 51’  P 52  ’  P 53’  P 61’  P 62  i j 

and  (11) 

dj  ={j,  u\u1,u2,v,v\  v„v2,w,  w\  w,  vv;,w2}^ 

The  Rayleigh-Ritz  approximation  yields  a  system  of  linear  equations  by  taking 
derivatives  of  the  variational  energy  equation  with  respect  to  the  independent  unknown  variables, 
as  follows: 


with 


a j(k) 


=  0 


(12) 


The  system  of  equations  is  then  expressed  in  the  matrix  form,  as  follows: 

(13) 


A 

C 

(*)  | 

[p] 

r_i 

Cr 

0 

1 

[dj 

t  1 

[f.J 

w 


A  =  J  Np  S  Np  dU 

U 

C  =  J  Np  BNd  dU 

U 

f1=0 

f2  =  J  Nj  o  dA 


(14) 


where  S,  B,  Np  and  Nd  are  matrices  for  the  compliance,  relationship  between  the  stresses  and 

displacements,  and  the  shape  functions  for  the  stress  and  displacement  degree  of  freedoms, 
respectively. 
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The  equations  for  the  constraint  conditions  at  the  interface  are  also  obtained  by 


— c—  =  0  and 

acf 


=o 


(15) 


which  yields  another  matrix  form,  as  follows: 


with 


1 

o 

■B 

o 

l 

O 

"S 

O 

_ i 

r  (kU 

P 

0  Qd  0  -Qd 

< 

d(t) 

-Qp  0  QP  0 

P(,) 

.  0  -Qd  0  Qd 

a 

(16) 


Qd  =«jNdr  N AdU 

U 

Qp  =tx2  {  Np  h  s  Np  dLJ 


(17) 


u 

where  h  is  a  matrix  containing  the  slopes  of  the  interfacial  surfaces  in  x-  and  y-directions  as  in 

equation  (4).  The  global  system  of  equations  is  then  formulated  by  combining  the  elemental 
stiffness  matrix  and  force  vectors  in  equation  (13)  and  equation  (16),  and  solved  numerically  to 
obtain  the  displacement  and  stress  results. 

1.2.2  Calculation  of  Effective  Elastic  Moduli 

Effective  elastic  moduli  of  the  woven  composites  are  calculated  by  solving  six 
different  cases  under  uniform  axial  and  shear  strain  loadings  ( ei ).  The  boundary  conditions  for 
each  loading  case  are  as  follows: 


(1)  For  e1=ex=Ua/Lx, 

u(0,  y,  z)  =0  ,  u{Lx,y,z)=U0 

v(x,0,  z)  =  v(x,  Ly,  z)  =  0  (18) 

w(x,  y,  0)  =  w(x,  y,  Lz)  =  0 
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(2)  For  F2  =£y  =V„/Ly  , 

u(0,  y,  z )  =  u(Lx,  y,  z)  =0 

v(x,0,z)=0  ,  v(x,  Ly,  z)  =  Va  (19) 

w(x,  y,  0)  =  w(x,  y,  L,  )=0 

(3)  For  e3=Et=W0/Lt, 

u(0,  y,  z)  =u{Lx,  y,  z)  =  0 

v(x,  0,  z)  =  v(x,  Ly  ,z)=0  (20) 

w(x,  y,0)  =0  ,  w(x,  y,Lz)=Wo 

(4)  For  =yyz  =VjLz  , 

u(0,  y,  z)  =u(Lx,  y,  z)  =  0 

v(x,y,0)  =  0  ,  v(x,y,Lz)  =  V0  (21) 

w(x,  y,  0)  =  w(x,  y,L_ )  =0 

(5)  For  e5  =YXZ  =  U„/LZ  , 

u(x,  y,0)=0  ,  u(x,  y,  L_ )  =  Ua 

v(x,  0,  z)  =  v(x,  Ly ,  z)  =  0  (22) 

w(x,  y,  0)  =  w(x,  y,Lz)  =  0 

(6)  For  =7^  =Ua/Ly  , 

u(x,  0,  z)  =0  ,  u(0,  Ly,  z)  =  U0 

v(x,  0,  z)  =  v(x,  Ly  ,z)=0  (23) 

w(x,  y,  0)  =  w(x,  y,Lz)=0 

where  U 0  ,Va  ,W0  are  the  uniformly  applied  displacement  components. 

For  each  uniform- strain  boundary  condition,  effective  stresses  (a,.)  are  calculated 
by  taking  a  volumetric  average,  as  follows: 


£ct,.(x,  y,z)  dxdydz 


C7  = 


y 


O'  =  i,  •  •  •  ,6) 


(24) 
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Components  of  6x6  effective  stiffness  matrix  ( [c (/  J )  and  effective  compliance  matrix  ( [st])  are 
then  obtained  by  the  following  equation: 


where 


W=[cjfe}  and  [sj=[cj  ' 


Cu 

r 

^ 12 

r 

^ 13 

r 

^ 14 

cl5 

CI6 

C21 

c 

22 

c 

^ 23 

c 

^ 24 

c25 

C26 

c 

^ 31 

c 

^  32 

c 

^33 

c 

^ 34 

c 

^ 35 

c 

^  36 

C4I 

c 

^ 42 

c 

^ 43 

c 

^44 

c45 

C46 

c5I 

c 

^ 52 

c 

^ 53 

c 

^  54 

C55 

C56 

C61 

c 

^ 62 

c 

^ 63 

c 

^ 64 

C65 

C66 

(25) 


(26) 


Finally,  the  three-dimensional  effective  elastic  moduli  are  calculated  by  the  following: 


E 


x 


G 


yz 


V 


xy 


(27) 


1.2.3  Calculation  of  Residual  Stresses 

When  the  laminated  composites  are  cured  and  cooled  to  room  temperature, 
residual  stresses  will  exist  because  thermal  contraction  of  each  ply  is  anisotropic.  When 
moisture  is  subsequently  absorbed,  hygro  expansion  is  also  anisotropic.  Therefore,  stress 
calculation  in  laminated  composites  should  include  both  the  mechanical  and  the  residual  stresses 
due  to  the  hygrothermal  effect.  The  variational  energy  term  for  the  hygro  thermal  effect  [1]  is 
calculated  as  follows: 


Fa)  = 


(28) 
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where  e-k)  =  {e\HgIobal  are  hygrothermal  strain  components  in  the  global  coordinate  system.  The 


hygrothermal  strain  components  ( e.k)  =  {e}on)  in  the  on-axis  coordinate  system  are  expressed  as 

e!k)  =  aA7  +  (3,c  (29) 

which  are  related  with  the  global  ones  by  the  following  equation: 

{<•)"„.,  =[rJoa]"r  [r,(9,)]r  {<■}",  (30) 

where  AT  is  a  temperature  difference  between  curing  and  operating  conditions,  and  c  is 

moisture  content  after  the  curing.  The  0,  is  the  angle  between  the  principal  material  direction 

and  the  yam  direction,  and  02  is  the  angle  between  the  yam  direction  and  the  global  coordinate 

direction.  The  [7,(0,)]  and  [7,  (02)]  are  tensor  transformation  matrices  for  fiber  orientation  and 

yarn-crimping  angles,  respectively,  as  follows: 

n  i  0  0  0  -2m1n1 

m]  0  0  0  2mInI 

0  10  0  0 
0  0  ml  11/  0 

0  0  —rij  m1  0 

m1n1  0  0  0  mf  —n] 

and  for  warp  yams,  the  following  applies: 

0  n\  0  2m2n2  0 

1  0  0  0  0 

0  m]  0  —  2  in ,  n ,  0 

0  0  m2  0  —n2 

0  m2n2  0  m22  —n]  0 

0  0  n2  0  m2 
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or  for  fill  yams,  the  following  applies: 

1  0  0  0  0 

0  m]  n\  —2m2n2  0 

0  n]  m  ]  2m2n2  0 

0  m2n2  —m2n2  m2  —  n  ;  0 

0  0  0  0  m2 

0  0  0  0  —n2 

where  m7  =  cost©,)  ,  n,  =sin(07)  ,  m2  =cos(02)  and  n2  =sin(02) . 

1.2.4  Calculation  of  Failure  Strength 

Failure  analysis  to  predict  the  critical  load  and  the  damage  mode  is  achieved  by 
applying  failure  criteria  to  the  reliably  calculated  stress  components,  not  by  interpolating  the 
displacement  results  as  in  the  displacement-based  method.  Assuming  that  all  stress  components 
increase  proportionally  to  the  applied  loading,  the  following  results: 

max  _  applied  (34) 

The  strength  ratio  (R)  can  be  split  into  mechanical  (  R"' )  and  residual  (R")  parts 
on  an  assumption  that  each  part  of  the  stresses  acts  independently,  so  that  equation  (34)  becomes 

a;,,ax  =R"’o”  +R’a;  (35) 

For  a  given  hygrothermal  combination  of  the  cure  temperature  and  the  moisture 
content,  the  residual  stresses  are  fixed.  When  mechanical  loads  are  applied  to  the  laminate,  the 
maximum  load  that  the  laminate  can  sustain  is  then  given  by  the  mechanical  part  of  the  strength 
ratio.  The  mechanical  strength  ratio  can  be  solved  by  letting  the  residual  strength  ratio  equal 
unity. 


[t2  (02 )] 


0 

0 

0 

0 

n  2 
m. 


(33) 
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1.2.4. 1 


Quadratic  failure  criteria 


Quadratic  failure  criteria  in  stress  space  consist  of  linear  and  quadratic 

invariants  as  follows: 


_ max _ max  .  n _ max  r  r\ 

pi 


(36) 


where  F..  and  F.  are  strength  parameters  in  stress  space.  By  letting  the  residual  strength  ratio 
equal  unity,  the  following  equation  results: 

am(Rm)2  +(bm  +  a )R +(ar  +br -1)  =  0  (37) 

where 


a  =  Fijai  CT;  ,  b  =Fiai 

amix  =  2Fijo'"orj 

ar  =  FiPp]  ,  br=Fp\ 


(38) 


In  this  study,  the  Tsai-Wu  failure  criterion  is  used  with  interaction  terms  of  F*2  =  Fn  =  -0.5  [6]. 

1.2. 4. 2  Maximum  stress  failure  criteria 

To  understand  the  basic  damage  mechanism  during  the  load  increase,  it 
is  important  to  identify  the  most  critical  stress  component.  To  distinguish  the  most  critical  stress 
component  from  the  others,  it  is  advantageous  to  use  maximum  stress  failure  criteria.  Note  that 
there  is  no  interaction  term  between  the  stress  components.  By  considering  both  the  mechanical 
and  the  residual  stresses,  the  failure  occurs  when  one  of  the  following  conditions  is  met: 

(1)  For  tensile  or  compressive  stresses  (i  =  l  ,2 ,3): 


if  a"'  +ar  >0  ,  Rm 

1  1  1  C5 m 

Xf  -a 


if  a'"  +o.  <0  ,  R;1  = 


1  1 


a 


(39) 
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(2)  For  shear  stresses  (i  =  4 ,5 ,6): 


R' 


s,-a; 


(40) 


where  Xf  ,  X-  and  S)  arc  tensile,  compressive  and  shear  strengths,  respectively. 


1.3  NUMERICAL  RESULTS  AND  DISCUSSION 

The  mixed  finite  element  method  is  implemented  into  an  in-house  computer  program, 
“3Dwoven.”  The  program  is  based  on  a  spreadsheet  with  user-friendly  input  and  output  routines. 

The  present  method  is  applied  to  the  analysis  of  flat  and  woven  laminated  composites. 
First,  displacements  and  stresses  of  these  composites  are  calculated  and  compared  with  analytic 
and/or  conventional  displacement-based  finite  element  solutions.  Second,  three-dimensional 
effective  elastic  moduli  are  calculated  for  several  flat/woven  laminated  composites,  and 
compared  with  existing  experimental/numerical  results.  Last,  first-ply  and  last-ply  failure  loads 
as  well  as  their  damage  modes  are  predicted  with  the  present  method.  The  numerical  predictions 
on  the  failure  strength  and  the  damage  mode  are  compared  with  experimental  results  that  were 
previously  observed  on  flat  laminated  composites  and  woven  model  laminates  with  one¬ 
dimensional  yarn  crimping. 

1.3.1  Flat  Laminated  Composites 

We  solved  a  class  of  boundary  value  problems,  known  as  the  free-edge  problem, 
in  which  a  flat  laminate  of  finite  width  is  subject  to  a  uniform  axial  displacement  ( Ua ).  The 
origin  of  coordinates  is  located  at  the  center  of  the  laminate,  and  the  laminate  is  symmetric 
(0(z)  =0(-z)).  Each  layer  is  treated  as  a  transversely  isotropic  material  with  a  lay-up  of  [0/90]s, 
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where  0°  is  parallel  with  the  x-axis,  as  in  Figure  2.  The  layers  are  of  equal  thickness,  h,  and  the 


laminate  width  is  2b  =16  h.  The  material  properties  are  listed  in  Table  1. 


Y 


Figure  2.  Flat  Laminated  Composites. 


Table  1 

Three-Dimensional  Properties  of  Unidirectional  T300/N5208  Composite 


Ex 

[GPa] 

Ey 

[GPa] 

Ez 

[GPa] 

Vvz 

GXy 

[GPa] 

Gxz 

[GPa] 

Gyz 

[GPa] 

181 

10.3 

10.3 

0.28 

0.28 

0.5 

7.17 

7.17 

7.05 

The  following  boundary  conditions  are  applied  to  simulate  a  tensile  loading 
subject  to  a  uniform  displacement  in  the  ^-direction.  The  axial  displacement  in  the  x-direction  at 
x=0  (yz-surface)  is  fixed,  and  at  x=Lx  is  prescribed  with  U0 .  The  symmetric  boundary  condition 
is  enforced  at  y=0  (xz- surface)  by  setting  v=0.  The  zero  vertical  displacement  at  z=0  simulates  a 
case  in  which  laminates  are  symmetrically  stacked. 


u\k)(0,  y,  z)  =  uf\0,  y,  z)  =  0 


u(k\0,y,z)=u*(k)  (0,  y,  z)  =  0 


u}k)(Lx,  y,  z)=u(2k\Lx,y,  z)=U 0  ,  ua)(Lx,  y,  z)  =  u*tk)  (Lx,  y,z)  =  ^y 

v\k)(x,  0,  z)  =v(k\x,  0,  z)  =  0  ,  v(k)(x,  0,z)  =v*ik)  (x,  0,  z)  =  0 
W(k)  (x,  y,0)=0 


(41) 
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The  penalty  parameters  in  equation  (6)  are  chosen  as  ay  =  10'  and  OC,  =  10 2  to 
avoid  numerical  instability.  Two  different  meshes  are  generated  for  the  present  mixed  method: 
having  the  same  number  of  divisions  in  the  y-axis,  one  has  only  one  sublayer,  and  the  other  has 
three  sublayers  in  each  ply  (subregion)  in  the  z-axis.  The  number  of  divisions  and  mesh  are 
shown  in  Figure  3  (a)  and  (b). 

Stress  and  displacement  results  are  compared  with  the  analytic  and  the 
displacement-based  finite  element  solutions.  The  analytic  solutions  are  obtained  by  a  two- 
dimensional  mixed  analysis  with  18  sublayers  suggested  by  Pagano  [3].  For  the  displacement- 
based  finite  element  method,  three-dimensional  eight-node  brick  elements  are  used  with  two 
different  meshes;  one  has  two  divisions  and  the  other  has  10  divisions  in  the  z-direction  in  each 
subregion,  as  in  Figure  3  (c)  and  (d).  The  interfacial  stresses  are  calculated  by  interpolating  the 
elemental  stresses  at  the  Gaussian  integration  points  into  the  nodal  points  along  the  interface. 
Thus,  two  normal  and  shear  stresses  are  calculated  by  interpolating  those  of  the  upper  and  the 
lower  elements  at  the  interface. 

Figure  4  shows  the  results  of  the  present  method  (mixed)  compared  with  the 
analytic  and  FEM  solutions.  The  results  show  that  the  normal  and  shear  stresses  become  singular 
at  the  free-edge  (  y  =  L  )  because  of  the  discontinuity  in  the  elastic  properties.  The  present  and 

analytic  methods,  which  are  both  the  mixed  methods,  yield  nearly  identical  results  for  the 
transverse  displacement  at  the  top  surface  [Figure  4  (a)],  the  normal  stress  along  the  [0/90] 
interface  [Figure  4  (b)],  and  the  normal  stress  along  the  central  surface  [Figure  4  (d)],  whereas  the 
displacement-based  FEM  shows  little  difference  with  them.  The  FEM  does  not  yield  an  accurate 
solution  without  a  sufficient  number  of  sublayers  in  the  z-direction,  whereas  the  present  method 
shows  an  excellent  agreement  even  with  one  layer  except  at  a  region  close  to  the  free  edge. 
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(a)  Mixed  ( Nx  =6,  Nv=24,  Nz=l ). 


(b)  Mixed  ( Nx  =4,  Ny  =24,  Nz=3 ). 


(c)  FEM  ( Nx  =  2  ,  Ny  =20 ,  Nz=2).  (d)FEM(A^v  =  7,  Ny  =  20,  Nz=10). 

Figure  3.  Geometry  and  Number  of  Divisions  for  Present  Mixed  and  Displacement-Based  Finite  Element 
Methods. 
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(a)  Transverse  displacement  at  top  surface  (  z  =  2h  ) 


4  - 


3  - 


— * — 

"Mixed  (Nz=l) 

— © — 

'  Mixed  (Nz=3) 

—  -A  — 

FEM  (Nz=2) 

-  -  O  - 

-FEM(Nz=10) 

Pagano  [3] 

1  - 


0  I WMI ItftK- Big  8 @(  BftlB gt  A - □  &■ 


0.2 


A 


£ 

A 

/  ti 

XJ 


0.4  0.6 

y/Lv 


0.8 


(b)  Distribution  of  a  along  [0/90]  interface  (  z  =  h ) 


Figure  4.  Stress  and  Displacement  Results  for  Flat  Laminated  Composites. 
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(c)  Distribution  of  xyz  along  [0/90]  interface  (  z  =  h ) 
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(d)  Distribution  of  cr  along  central  plane  (z  =  0) 


Figure  4.  Stress  and  Displacement  Results  for  Flat  Laminated  Composites  (concluded). 
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While  the  analytic  solution  yields  zero  shear  stress  with  a  high  peak  at  the  free 


edge  (y  =Ly ),  the  present  and  FEM  solutions  give  finite  values,  as  Figure  4  (c)  shows.  This  is 
because  the  linear  shape  functions  used  in  both  finite  element  methods  are  not  accurate  enough 
to  capture  the  drastic  stress  change  at  the  free  edge.  Although  the  increment  of  the  sublayer  in 
the  present  mixed  method  makes  the  peak  value  higher,  it  creates  wiggles  in  the  shear-stress 
distribution  near  the  free  edge.  This  is  because  the  high  stress  gradient  at  the  free  edge 
influences  the  stress  field  inside  the  edge.  As  pointed  out  earlier,  the  excessive  continuity  for 
stresses  should  be  avoided  at  singularities  and  at  abrupt  material  property  change  interfaces. 
Therefore,  the  penalty  method  is  more  suitable  than  the  irreducible  formulation  because  it  can 
relieve  the  excessiveness  by  controlling  the  penalty  parameter  for  the  stress  constraint  condition 
(oc2 ).  Note  that  in  this  case  of  such  an  extremely  high  stress  gradient,  even  the  penalty  method 
cannot  cure  the  problem  completely. 

The  normal  stress  at  the  central  surface  [Figure  4  (d)]  with  one  sublayer  shows 
good  agreement  with  the  analytic  solution  except  for  a  hump  at  the  free  edge.  This  hump  does 
not  appear  with  three  sublayer  solutions. 

1.3.2  RVE  of  Woven  Fabric  Composites 

The  RVE  of  the  model  is  divided  into  several  subregions;  each  subregion  is 
occupied  by  a  characteristic  fabric  yarn  or  a  matrix  (see  Figure  5).  The  L  and  Lv  variables 

represent  the  length  of  RVE  in  the  x-  (warp)  and  y-  (fill)  directions,  and  tw  and  t  f  half  of  the 

thickness  of  the  warp  and  fill  yarn,  respectively.  The  yarn  is  assumed  as  transversely  isotropic, 
and  the  matrix  as  isotropic  materials.  Each  yarn  and  the  matrix  subregion  of  the  RVE  are 
discretized  into  several  finite  elements  in  the  longitudinal  and  transverse  directions.  The  Nx  and 
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Figure  5.  Representative  Volume  Element  of  a  Plain-Weave  Composite.  Numbers  in  circles  indicate  the 
numbers  of  subregions. 

Ny  variables  represent  the  numbers  of  subdivisions  in  half  of  the  length  of  RVE  in  the  x-  and  the 
y-  directions  ( Lx  /2  and  Ly  j2 ),  respectively. 

The  cross-sectional  boundary  of  the  yam  is  confined  by  hL  =  ht  (lower  boundary) 
and  hu  =  h2  (upper  boundary).  Because  of  yam  waviness  and  the  elliptical  cross-sectional 

boundary  of  the  yarns,  hL  and  hu  are  functions  of  both  x  and  y  .  Yam  waviness  at  each 
subregion  is  assumed  to  be  sinusoidal  functions.  Lower  and  upper  surface  coordinates  of  the 
yarn  subregions  are  as  follows: 
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where  a  superscript  indicates  the  subregion  number.  Lower  and  upper  surface  coordinates  of  the 
bottom  and  top  matrix  subregions  are  as  follows: 

(1)  for  bottom  matrix  subregion  (subregion  5), 

h‘i"  =  0,  and 


h™  = 


t j (1  —  cos——)  +  tw(l  —  cos——)  for  0<x<Lx/2,0<y<Ly/2 
Lx  Ly 

t  f  (1  +  COS  ^-)  +  tw  (1  -  cos  J-)  for  Lj  2<x<  Lx  ,  0  <y<  Ly  /2 
t  f  (1  -  cos  — )  +  tw  (1  +  cos  — )  for  0  <  x<  Lx/2  ,  Ly/2<  y  <  Ly 


(43) 


^(1  +  cos— )  +  fw(l  +  cos— )  for  Lx/2<x  <Lx  ,  Ly/2<y  <Ly 

I j  Z. 


(2)  for  top  matrix  subregion  (subregion  6), 
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(44) 


The  following  boundary  conditions  are  prescribed  to  simulate  a  tensile  loading 
subject  to  a  uniform  displacement  in  the  .r-direction  with  lateral  constraint  in  the  y-direction.  The 
zero  vertical  displacement  at  z=0  simulates  a  case  in  which  two  RVEs  are  symmetrically 
stacked,  as  shown  in  the  following  equation: 
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uf  \0,  y,  z)  =u(2k)(0,  y,  z)  =  0  ,  u(k> (0,  y,  z)  =  u*{k)  (0,  y,  z)  =0 


<\Lx,  y,  z)  =  u[k)(Lx ,  y,  z)=U0  ,  u(k)(Lx,  y,z)  =  u  *<*>  (L,  y,z)  =  ^ 

vf\x,0,  z)  =  v<2k)(x,  0,  z)  =  vi1k)(x,  Ly,z )  =vf)(x,  Ly,  z)  =  0  (45') 

v(k> (x,  0,  z)  =  v*(k)  (x,0,  z)  =  vik)(x,  Ly,  z)  =  v*ik)  (x,  Ly,z )  =0 
w<k)(x,y,  ())=() 

Figure  6  (a)  shows  a  deformed  shape  under  the  above  boundary  conditions.  The 
top  surface  of  the  RVE  is  twisted  because  of  its  antisymmetric  geometry  in  the  x-  and  y- 
directions.  Figure  6  (b)  shows  the  antisymmetric  distributions  of  the  vertical  displacement  at  the 
intersection  of  the  top  surface  and  the  xz-planes  at  y  =  0  and  y  =  Ly . 

The  thickness  of  the  matrix  subregions  (subregions  5  and  6)  at  four  comer  points 
is  zero  according  to  the  model.  Physically,  the  lower  (  w, )  and  upper  ( w2 )  vertical 
displacements  at  these  corners  should  be  the  same.  However,  because  of  the  numerical  errors, 
they  do  not  match  with  each  other  with  the  coarse  meshes  ( iVv  <  3 ),  as  in  Figure  6  (c). 
Therefore,  finer  meshes  ( Nx  >4)  should  be  used  to  achieve  the  interfacial  continuity,  and 
N  x  =  6  is  chosen  in  this  study. 

While  one  displacement  penalty  parameter  is  set  as  CL,  =  103 ,  two  stress  penalty 

parameters  are  chosen  as  a2  =  102  and  a,  =0  for  a  sensitivity  study.  The  latter  case  (a,  = 0 ) 
means  no  stress  constraint  condition  is  enforced.  Figure  6  (d)  shows  that  the  vertical 
displacement  distributions  are  almost  identical  with  two  different  a2 ,  which  indicates  that  the 
stress  continuity  condition  has  a  negligible  influence  on  the  displacement  results. 
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(a)  Deformed  shape  ( Nx  =6) 
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(b)  Vertical  displacements  at  top  surface 


x/Lx 

0  0.2  0.4  0.6  0.8  1 


(d)  Vertical  displacements  with  two  a2 


Figure  6.  Displacement  Results  of  RVE  of  Woven  Composites. 


Figure  7  shows  the  normal  and  shear  stress  distributions  along  an  interfacial  line 
in  Figure  5.  These  interfacial  stresses  are  the  local  ones  described  in  equation  (4),  which  are 
transformed  from  the  stress  components  in  the  global  coordinate  system  by  the  slopes  of  the 
interfacial  surfaces.  The  subscript  ( k )  indicates  the  bottom  matrix  subregion  (subregion  5),  and 
(/)  indicates  the  upper  yams  lying  on  top  of  the  matrix  (i.e.,  subregion  1  at  0  <  x<  Lj2  and 
subregion  3  at  Lx  /2  <  x  <  Lx ).  Figure  7  (a)  and  (b)  show  that  the  interfacial  normal  stresses 
from  the  lower  and  the  upper  subregion  agree  well  with  each  other,  with  only  one  sublayer  in  the 
thickness  (z)  direction.  The  normal  stress  continuity  can  be  achieved  well  even  without  the  stress 


26 


x/Lx 


(c)  Shear  stress  ( xxz )  with  a,  =100 


(d)  Shear  stress  (xw )  with  a2  =  0 


Figure  7.  Interfacial  Normal  and  Shear  Stress  Distributions  of  RVE  of  Woven  Composites  with  Two 
Different  Penalty  Parameters  for  Stress  Continuity  Condition. 


constraint  condition  ( oc,  =  0 ).  It  also  shows  a  smooth  transition  of  the  stress  distribution  with  a 
significant  change  in  the  material  properties  at  x  =  Lj2  . 

Figure  7  (c)  and  (d)  show  that  the  interfacial  shear  stress  continuity  is  achieved 
fairly  well  with  the  present  method,  except  the  region  near  x  =Lj2  ,  where  the  high  stress 
gradient  is  observed.  The  reason  for  the  high  stress  gradient  in  the  local  shear  stresses  is  that  the 
local  shear  stresses  ( <34  and  d5 )  are  highly  affected  by  the  global  axial  stresses  (  ct;  and  ct2  )  as 
indicated  in  equation  (4),  and  these  axial  stresses  change  abruptly  with  the  change  in  the  material 
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properties  at  this  region.  The  interfacial  shear  stresses  do  not  match  well  at  this  region  because 
the  thickness  of  subregion  3  is  zero  at  x  =  Lj2  .  While  two  subregions  (subregions  5  and  1)  are 

considered  in  calculating  the  interfacial  stresses  at  the  left  side  of  x  =  Lj2  ,  three  subregions 
(subregions  5,  3  and  1)  are  considered  in  the  calculation  at  the  right  side  because  of  the  zero 
thickness  of  subregion  3.  Therefore,  the  stress  continuity  condition  becomes, 

(1)  at  left  side  of  *  =  Lx/2  , 

P£=P$  and  P™  =  P™  (46) 


(2)  at  right  side  of  x  =  Lx/2  , 

p<5>  _  p(3)  _  p(,3)  _  pt?) 
r42  r  41  r42  r  41 


and 


p(5)  _  pi.  3)  _  p  (3) 
r  52  r  51  r  52 


p(D 
r  51 


(47) 


where  Pt(k>  is  the  local  component  of  the  interfacial  stresses  at  the  k-th  subregion.  However,  it  is 

hard  to  satisfy  such  a  continuity  condition  with  the  zero  thickness  because  of  the  numerical  error 
in  evaluating  the  stress  components.  The  numerical  error  in  the  axial  stresses,  whose  magnitudes 
are  much  larger  than  those  of  the  shear  stresses,  affects  the  interfacial  shear  stresses  significantly, 
so  that  jumps  and  mismatches  are  observed  at  this  region. 

Figure  7  (c)  and  (d)  also  show  that  the  shear  stress  distribution  is  smoother 
without  the  stress  constraint  condition  ( oc2  =  0 )  than  a 2  =  100 .  As  observed  in  the  flat- 
laminated  case,  the  excessive  stress  continuity  conditions  are  not  necessary  in  the  present  mixed 
method,  and  should  be  avoided  at  the  stress  singularity  or  the  material  mismatch.  Not  shown  on 
the  figure  are  the  results  for  cl2  » oc7 ,  which  make  a  little  improvement  in  the  stress  continuity 
but  cause  the  displacement  results  to  be  unrealistic  and  far  different  from  the  one  in  Figure  6  (a). 
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1.3.3  Effective  Elastic  Moduli  of  Plain-Weave  Laminates 


Zhang  and  Harding  [7]  used  a  strain  energy  equivalency  principle  to  predict  the 
elastic  properties  of  a  plain-weave  composite.  The  finite  element  method  was  used  to  evaluate 
the  strain  energies.  They  applied  this  method  to  a  one-directional  undulation  model  in  the 
loading  direction  only.  Comparisons  were  made  with  experimental  data  for  a  plain-weave 
carbon  epoxy  laminate  [8].  Because  of  the  one-dimensional  model,  discrepancies  occurred  for 

the  in-plane  shear  modulus  ( Gn  )  and  the  properties  in  the  transverse  direction  (  E  and  v ). 

Naik  and  Ganesh  [9]  suggested  two  refined  models,  slice  array  model  (SAM)  and 
element  array  model  (EAM),  and  also  suggested  modifications  of  the  existing  simple  models, 
modified  mosaic  parallel  model  (MMPM)  and  modified  Kabelka’s  model  (MKM).  These 
models  predicted  two-dimensional  elastic  properties,  considering  the  actual  yarn  cross-section 
geometry,  possible  gap  between  two  adjacent  yarns,  and  undulation  and  continuity  of  yams 
along  both  warp  and  fill  directions.  The  effective  moduli  are  calculated  by  the  various  models 
for  plain- weave  composites  with  E-glass/epoxy  and  carbon/epoxy  materials. 

Figure  8  (a)  shows  an  RVE  of  the  present  model  with  the  maximum  yam 

thickness  (a)  and  the  wavelength  of  the  yam  (A.),  whose  ratio  (a  IX)  represents  the  waviness 
ratio.  The  overall  volume  fraction  filled  with  yam  subregions  in  the  RVE  is  approximately  0.64. 
The  moduli  agree  fairly  well  with  both  the  experimental  and  the  numerical  results  for  various 
waviness  ratios,  as  shown  in  Figure  8.  Not  plotted  in  the  figure,  the  present  method  can  calculate 

three-dimensional  moduli  and  Poisson’s  ratios,  such  as  E_ ,  G  _ ,  G  .  V~  and  v"  ,  as  opposed  to 
the  existing  two-dimensional  methods. 
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(a)  RYE  and  dimension 


(b)  Ishikawa  [8]  (a /  A=0.025) 
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(c)  Naik  [9]:  E-glass/epoxy  (a/ X  =0.075) 


(d)  Naik  [9] :  carbon/epoxy  (aA,=0.035) 


Figure  8.  Comparison  of  Numerical  Prediction  of  Effective  Elastic  Moduli  with  Existing  Results. 


1.3.4  Failure  Analysis  of  Model  Laminated  Composites 

In  situ  damage  observation  is  made  with  flat  and  model  laminates  containing  one- 
directional  yam  crimping  [10].  Test  specimens  were  loaded  in  tension  in  the  x-direction  (warp 
direction)  in  a  portable  load  frame  placed  on  a  microscopic  stage.  Laminates  were  made  from 
AS4/3501-6  graphite/epoxy  unidirectional  prepreg,  whose  properties  are  listed  in  Table  2  [11]. 
The  present  numerical  prediction  is  compared  with  the  experimental  results  for  a  cross-ply  flat 
laminate,  [90/0]2S,  and  a  model  laminate,  ([9O2/O2L  ,  [O2/9O2L  ,  0.050).  The  notation  used  for  the 
model  laminate  indicates  the  lamination  sequence  away  from  the  wavy  region,  midsection 
lamination  sequence  of  the  wavy  region  and  the  waviness  ratio,  respectively.  The  hygrothermal 
conditions  of  AT  =  -95  °C  and  c=0.005  are  used  in  the  calculation. 
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Table  2 

Three-Dimensional  Properties  of  Unidirectional  AS4/3501-6  Composite 


Engineering  Constants 

Ex 

Ey 

Ez 

v« 

V.v, 

Gyz 

Gxz 

Gxy 

[GPa] 

[GPa] 

138 

10 

10 

0.3 

0.3 

0.53 

2.9 

5.5 

5.5 

Strength  Data  [MPa] 

Tension 

Compression 

Shear 

XTx 

XT, 

X\ 

xcx 

Xy 

K 

s* 

1930 

52 

52 

1450 

210 

258 

103 

93 

93 

In  both  the  flat  and  woven  cases,  the  first-ply-failure  (FPF)  occurs  at  the  90°  fill  yams  in  the 
form  of  transverse  matrix  cracking.  In  the  flat  laminate,  the  maximum  stress  occurs  at  y  =  Ly 

because  of  the  singularity  in  ov  at  the  free  edge.  The  stress  distribution  in  the  /-direction  is 
calculated  with  the  present  method  and  compared  with  an  analytic  solution  obtained  by  a  method 
suggested  by  Pagano  [3],  as  Figure  9  shows.  The  figure  clearly  shows  the  singularity  at  the  free 
edge.  Similar  distribution  is  observed  in  the  woven  model  laminates.  Therefore,  fine  meshes  are 
required  near  the  free  edge  for  an  accurate  failure  prediction. 


Figure  9.  6 x  Distribution  in  /-Direction  at  an  Interface  between  Lower  0°  and  Lower  90°  Ply. 
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In  the  model  laminates,  the  calculated  stress  distribution  in  the  v-direction  shows 


that  the  maximum  stress  occurs  at  the  upper  fill-yam  near  the  crimping  at  x  =  0.57  L  .  The 
calculated  location  agrees  well  with  the  in  situ  observation  of  the  first  matrix  cracking,  as  Figure 
10  shows. 


First  crack  at  y=0.54  Ly 


Figure  10.  Cross  Section  of  Woven  Model  Laminate  when  the  First  Cracking  is  Observed  [10]. 

A  primitive  progressive  damage  model  is  achieved  by  degrading  transverse  yam 
modulus  with  a  matrix  degradation  factor  of  E*m  =  0.2  .  The  last-ply-failure  (LPF)  occurs  at  the 

0°  warp  yarns  in  the  form  of  fiber  breakage.  The  strengths  at  FPF  and  LPF  are  plotted  in  Figure 
11  with  fairly  good  agreement  with  the  experiment.  The  Tsai-Wu  quadratic  failure  criterion 
provides  more  conservative  and  better  results  than  the  maximum  stress  criterion  because  the 
former  considers  interactions  between  the  stress  components.  However,  the  maximum  stress 
criterion  is  useful  to  identify  the  critical  stress  component  for  the  damage  mode  prediction. 

1.4  SUMMARY 

Three-dimensional  displacements  and  stresses  are  analyzed  numerically  based  on 
Reissner’s  mixed  variational  principle.  The  three-dimensional  model  is  treated  semi-two- 
dimensionally  by  making  an  assumption  on  the  interlaminar  stress  variations  and  integrating  the 
variational  energy  in  the  thickness  direction.  Additional  energy  terms  are  added  to  the 
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(a)  Flat:  ([90/0]2s) 


2  -i 


FPF  LPF 


(a)  Woven:  ([902/02]s ,  [02/902]s ,  0.050) 


Figure  11.  Numerical  Prediction  of  Strains  at  FPF  and  LPF.  Total  strain  is  a  sum  of  mechanical  and 
residual  strains. 


variational  energy  to  impose  the  displacement  and  stress  continuity  at  the  interface  by  the  penalty 
approach.  Two  penalty  parameters  are  employed  to  enforce  the  displacement  and  the  stress 
continuity  conditions,  respectively.  The  Rayleigh-Ritz  approximation  with  polynomial  shape 
functions  yields  a  system  of  linear  equations  by  taking  derivatives  of  the  variational  energy 
equation  with  respect  to  the  independent  unknown  variables. 

The  present  method  is  applied  to  analyze  flat  laminated  composites  with  a  free  edge  and 
the  RVE  of  the  plain- weave  composites.  The  results  are  compared  and  validated  with  the 
displacement-based  FEA  and/or  analytic  solution.  Since  the  stresses  are  evaluated  pointwise 
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without  any  interpolation  of  the  displacement  results,  more  accurate  interlaminar  stresses  are 
obtained  at  the  interfaces  between  two  different  materials  with  few  sublayers  compared  with  the 
displacement-based  FEA. 

The  interfacial  normal  and  shear  stress  continuity  is  achieved  well  with  the  penalty 
approach,  except  in  the  region  in  which  the  thickness  of  the  subregions  is  small.  It  is  found  that 
the  imposition  of  the  displacement  continuity  condition  is  more  important  than  that  of  the  stress 
continuity  condition.  Furthermore,  the  excessive  continuity  condition  in  the  stress  fields  is  not 
necessary  and  may  induce  convergence  instability  in  the  data  of  the  displacement  as  well  as  the 
stress  fields.  Imposing  only  the  displacement  constraint  without  the  stress  constraint  yields  a 
smoother  interfacial  normal  and  shear  stress  distribution  than  the  case  that  considers  both 
constraint  conditions. 

The  reliable  stress  calculation  is  used  to  predict  the  effective  elastic  moduli.  The 
numerical  calculation  with  the  present  model  of  the  RVE  shows  good  agreement  with  previous 
experimental  and  numerical  results  made  for  both  flat  and  woven  laminated  composites  with 
various  yarn-waviness  ratios.  Since  the  present  method  calculates  the  three-dimensional  elastic 
moduli  with  three-dimensional  geometry,  it  can  be  used  not  only  for  the  2-D  but  also  for  3-D 
woven  composites. 

Discrete  damage  analysis  is  achieved  with  the  reliable  prediction  of  the  stress  field  with 
the  existing  failure  criteria.  The  failure  analysis  includes  residual  stress  calculation  to  consider 
hygro thermal  effects.  The  damage  analysis  results  in  a  good  prediction  of  the  magnitude  and 
location  of  the  failure.  A  primitive  progressive  damage  analysis  is  established  by  degrading  the 
material  properties  of  the  damaged  yams  to  predict  the  final  failure  beyond  the  FPF. 
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With  the  presence  of  stress  singularity  at  the  free  edge,  the  numerical  calculation 
becomes  mesh-dependent.  Therefore,  it  is  desirable  to  eliminate  the  singularity  or  reduce  the 
dominance  of  the  singular  stresses  in  order  to  study  the  effect  of  the  yarn  crimping  on  the  failure 
of  the  woven  composites.  Non-straight-edge  specimens,  such  as  a  cruciform,  can  be  used  for  this 
purpose  to  make  the  most  of  the  stresses  carried  not  by  the  edges  but  by  the  middle  of  the 
specimen. 

1.5  CONCLUSIONS  AND  RECOMMENDATIONS 

We  have  developed  a  three-dimensional  model  for  the  three-dimensional  stress  analysis 
of  woven  fabric  composites.  The  model  yields  an  accurate  three-dimensional  displacement  and 
stress  solution  of  woven  fabric  composites  under  any  of  the  in-plane  and  the  out-of-plane  loading 
conditions  including,  but  not  limited  to,  extension,  bending  and  twisting.  The  model  can  obtain 
the  three-dimensional  effective  stiffness  matrix  for  woven  composites  that  a  designer  can  plug 
into  for  finite  element  structural  analysis.  The  model  can  also  make  an  accurate  damage 
prediction  for  woven  composites  through  accurate  in-plane  and  interlaminar  stress  calculations. 

This  present  mixed  formulation  offers  better  accuracy  at  high  stress  gradients  over  the 
current  state-of-the-art  displacement-based  method.  However,  by  introducing  more  variables, 
this  new  method  becomes  resource-intensive  in  terms  of  computer  memory  and  computational 
time.  A  future  effort  will  be  made  to  develop  a  model  based  on  the  displacement-based  method 
to  make  the  computation  economical  and  efficient. 
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2.  ANALYTICAL  CHARACTERIZATION  OF  GRAPHITIC  FOAMS 


2.1  INTRODUCTION 

In  recent  years,  there  have  been  an  increasing  number  of  applications  requiring 
lightweight  and  more  efficient  thermal  management,  such  as  high-density  electronics,  hybrid 
diesel-electric  vehicles,  communication  satellites,  and  advanced  aircraft.  The  primary  concerns 
in  these  thermal  management  applications  are  high  thermal  conductivity,  low  weight,  low 
coefficient  of  thermal  expansion,  high  specific  strength,  and  low  cost  [12].  Carbon  foam  was 
shown  to  demonstrate  numerous  unique  properties  that  make  it  an  attractive  material  for  use  of 
low-cost,  lightweight,  insulating,  energy-absorbing  structural  components.  Unique  properties  of 
the  carbon  foam  material  include  [13]: 

(1)  Precursors:  coal  extracts  are  inexpensive  (less  than  $0.10  per  kg)  and  readily  available. 

(2)  Manufacturing  of  the  foam  can  be  readily  scaled  up  by  continuous  extrusion  of  constant 
cross-section  parts,  or  net-shape  batch  production  of  special  shapes.  Required 
manufacturing  equipment  is  commercially  available;  projected  finished  material  cost  is 
less  than  $14  per  kg. 

(3)  Low  bulk  thermal  conductivity:  less  than  1.0  W/mK,  but  potential  for  high  thermal 
conductivity  if  the  foam  is  converted  to  graphite  through  heat  treatment  at  >2000°C. 

(4)  Fire  resistance:  once  carbonized  at  >1000°C,  the  foam  does  not  contain  a  sufficient 
volatile  material  with  which  to  support  combustion. 

(5)  It  will  not  give  off  noxious  or  hazardous  fumes  when  heated. 

(6)  Its  properties  can  be  readily  engineered  to  meet  different  requirements.  By  varying  the 
processing  conditions,  the  density,  compressive  strength,  and  ability  to  absorb  energy  can 
be  tailored  to  meet  specific  requirements. 
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(7)  Integration  with  other  materials:  examples  include  impregnation  with  phenolic  or  other 
resins,  lamination  with  Kevlar  tape,  and  lamination  with  a  phenolic -resin  skin.  Attaching 
fiber-reinforced  polymer  or  metallic  facesheets  allows  joining  to  other  components  by 
more  conventional  methods,  protects  the  foam  from  localizing  damage  or  abrasive  wear, 
and  transfers  loads  uniformly  to  the  foam. 

(8)  Machinability:  easily  cut,  milled,  turned,  etc.,  with  conventional  equipment  and  tooling. 

(9)  Formability:  foam  assumes  shape  of  mold  in  batch  operation  and  may  be  continuously 
formed. 

(10)  Joining:  using  a  coke  fusion  process.  This  feature  enables  foams  with  different 
mechanical,  thermal  or  electrical  conductivity  properties  to  be  joined  to  produce  a  highly 
tailorable,  anisotropic,  sandwich  material,  as  well  as  to  allow  repair  to  damaged 
structures. 

(11)  Impact  absorption:  carbon  foam  performs  better  than  conventional  polyurethane  foams 
that  are  currently  used  extensively  for  impact  absorption  in  aircraft. 

(12)  Additional  improvement:  additives  such  as  chopped  fibers,  nanofibers  or  nanotubes,  and 
crushed  calcined  cokes  can  add  significantly  to  the  strength  and  tailorability  of  the  foams; 
unidirectional  expansion  of  the  foam  and  the  orientation  of  fibers  within  the  matrix 
enable  the  production  of  anisotropic  foams  with  directional  properties. 

The  carbon  foam  macroscopically  possesses  an  isotropic  material  property.  However,  a 
microstructure  of  an  open-cell  foam  possesses  a  pentagonal  dodecahedron  structure  of  the  foam 
ligaments  oriented  approximately  109.47°  with  each  other.  The  microstructure  of  foams  reflects 
their  method  of  preparation,  which  usually  involves  a  continuous  liquid  phase  that  eventually 
solidifies.  Surface  tension  and  the  elementary  features  of  the  liquid  foam  structure,  which  are 
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required  to  minimize  surface  energy  during  the  foaming  process  (i.e.,  bubble  nucleation  process), 
results  in  three  films  that  always  meet  at  equal  angles  of  120°  to  form  a  film  junction  region 
called  a  plateau  border,  and  four  plateau  borders  always  join  at  the  tetrahedral  angle  of 

cos1  (1/ 3)  ~  109.47°  [14,15], 

Because  of  the  tetrahedral  cell  microstructure,  the  macroscopic  properties,  such  as  foam 
moduli  and  strengths,  are  critically  influenced  by  the  deformation  characteristics  of  the  cell 
ligaments.  Therefore,  to  develop  a  basic  understanding  of  the  performance  of  open-cell  foam 
materials,  it  is  critical  that  the  deformation  and  failure  mechanism  of  the  cell  ligaments  critically 
be  studied. 

Preliminary  research  reveals  that  the  graphitic  alignment  of  the  cell  ligaments  varies 
along  its  longitudinal  direction.  Processing  parameters,  such  as  temperature,  pressure,  etc., 
determine  the  porosity  and  the  graphitic  alignment  of  the  carbon  foam,  which  in  turn  determine 
its  geometries  and  material  properties.  Once  a  research  effort  investigates  appropriate 
microstructural  characterization  techniques  to  correlate  the  foam  microstructure  with  the 
processing  parameters,  the  foam  microstructural  geometry  and  the  material  properties,  including 
mechanical  elastic  moduli,  Poisson’s  ratio  and  thermal  conductivity,  etc.,  will  be  used  for  the 
mechanical  and  thermal  analysis. 

Because  of  slenderness,  each  ligament  can  be  considered  as  a  beam,  and  the  tetrahedral 
cell  microstructure  with  four  ligaments  as  a  frame  structure.  The  four  beams  are  located  in  three- 
dimensional  space.  The  cross  section  of  the  beam  varies  in  size  and  material  properties  in  the 
longitudinal  direction  along  the  ligament.  Because  of  the  complex  geometry  and  anisotropic 
material  properties,  it  is  appropriate  to  perform  the  analysis  numerically  to  obtain  accurate 
displacement  and  stress  field  solutions  of  foams.  The  numerical  analysis  will  predict 
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longitudinal  and  transverse  displacements  as  well  as  rotations,  and  calculate  the  reliable  stress 
and  strain  distributions  along  the  beam  ligaments  that  are  connected  and  located  in  the  three- 
dimensional  space  and  are  deformed  under  arbitrary  loading  conditions. 

2.2  PRELIMINARY  ANALYSIS  OL  LOAM  MODEL 
2.2.1  Generation  of  Unit  Cell  of  Carbon  Loam 

Because  of  the  randomness  and  complexity  of  the  microstructure  of  the  carbon 
foam,  it  is  difficult  to  consider  every  cell  of  the  whole  foam  individually.  Instead,  an  assumption 
is  made  that  one  of  the  cells  can  represent  a  certain  behavior  of  the  whole  foam  structure.  The 
microstructural  characterization  of  the  representative  cell  ligaments  will  then  be  correlated  with 
the  macroscopic  bulk  properties  by  a  statistical  approach.  A  series  of  databases  will  be  collected 
for  various  size  and  spatial  orientation  of  the  cell  ligaments,  as  well  as  for  the  property  variation 
due  to  the  graphic  alignment  along  the  longitudinal  direction  of  the  ligaments. 

The  first  step  in  the  structural  analysis  of  the  foam  is  to  generate  the  geometry  of 
the  unit  cell  of  the  foam  systematically.  The  unit  cell  can  be  generated  by  manipulating 
geometric  entities,  such  as  keypoints,  lines,  areas  and  volumes.  The  ANSYS  package,  a 
commercially  available  finite  element  package,  is  used  to  manipulate  those  entities  by  adding, 
subtracting  or  merging  methods. 

The  first  step  is  to  create  a  cube  with  dimensions  of  2ax2ax2a  in  x-,  y-  and  z- 
directions,  as  shown  in  Figure  12.  The  origin  (point  1)  is  located  in  the  center  of  the  cube.  The 
second  step  is  to  select  four  comer  points  (e.g.,  points  3,  5,  6,  8  in  Figure  12),  which  are  located 
diagonally  with  each  other  on  the  faces  of  the  cube.  Connecting  the  four  corner  points  then 
generates  a  tetrahedron,  whose  volume  is  VTetra  =  8/ 3 -a3  .  This  tetrahedron  can  be  divided  into 
four  subtetrahedra  that  contain  the  origin  (1-3-8-5,  1-8-6-5,  1-5-6-3  and  1-3-8-6).  Local 
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Figure  12.  Keypoints  and  Lines  for  Generating  a  Unit  Cell  of  Carbon  Foam, 
coordinate  systems,  whose  ^-directions  are  parallel  to  the  longitudinal  directions  of  the 
ligaments,  are  defined  by  using  four  lines  that  connect  the  keypoints  (1-3,  1-5,  1-6  and  1-8).  The 
local  coordinate  systems  are  useful  because  the  anisotropic  material  properties,  with 
consideration  of  the  graphitic  alignment,  are  input  easily  by  defining  the  longitudinal  directions 
of  the  ligaments. 

The  next  step  is  to  generate  four  spheres  on  the  four  corner  points  of  the 
tetrahedron.  The  spheres  represent  bubbles  that  are  produced  during  the  foaming  process.  The 
radii  of  the  spheres  will  determine  the  porosity  of  the  unit  cell  of  the  foam.  Figure  13  shows  the 
tetrahedron  and  the  spheres. 

The  next  step  is  to  subtract  the  spheres  from  the  tetrahedron  by  geometric 
manipulation.  The  remaining  media  is  the  unit  cell  of  the  foam.  The  volume  of  the  unit  cell 
VCM  is  calculated  automatically  by  the  ANSYS.  The  porosity  ((f))  of  the  foam  is  thus  calculated 
by  <f>  =  1  -VCeU/VTetra  .  Figure  14  shows  the  unit  cells  with  various  porosities,  and  Figure  15 
shows  one  of  the  ligaments  of  the  unit  cell  of  the  foam  at  different  view  angles. 
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Figure  13.  Tetrahedron  and  Spheres  to  Generate  a  Unit  Cell  of  Carbon  Foam. 


(a)  <|>  =  10%  (b)  <|>  =  78%  (c)  <|>  =  98% 

Figure  14.  Unit  Cells  of  Carbon  Foam  with  Various  Porosities. 


For  numerical  analysis,  it  is  advantageous  to  partition  the  ligaments  along  the  longitudinal  direc¬ 
tions  because  of  the  varying  material  properties  due  to  the  graphitic  alignment.  It  is  then  assumed 
that  material  properties  change  segment  by  segment.  The  number  of  partitions  will  determine  the 
accuracy  of  the  analysis.  The  ligament  partitioning  can  be  easily  generated  by  manipulating  the 
tetrahedra  at  the  beginning  step.  The  tetrahedra  can  be  duplicated  and  contracted,  leaving  the  origin 
in  the  same  location.  The  unit  cell  with  the  partitioned  ligaments  is  then  generated  by  following  the 
same  steps  as  generating  the  spheres  and  subtracting  them  from  the  partitioned  tetrahedra 
geometrically,  as  Figure  16  shows.  By  setting  mesh  parameters  for  a  desired  refinement,  and 
running  an  automatic  mesh  routine  of  ANSYS,  three-dimensional  tetrahedral  finite  elements  are 
automatically  generated  as  shown  in  Figure  17. 
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Figure  15.  A  Ligament  of  a  Unit  Cell  of  Carbon  Foam  at  Different  View  Angles. 


Carbon  Foam 


Figure  16.  A  Unit  Cell  of  Carbon  Foam  Partitioned  for  Varying  Material  Properties  Along  Ligaments. 
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Figure  17.  Finite  Element  Meshes  for  Unit  Cell  of  Carbon  Foam. 

2.2.2  FEA 

For  the  FEA,  appropriate  loading  and  boundary  conditions  should  be  specified  on 
the  unit  cell  of  the  foam.  Under  a  certain  loading  condition,  the  carbon  foam  behaves  according 
to  overall  material  properties,  such  as  bulk  modulus  and  Poisson’s  ratio.  The  overall  material 
properties  can  be  assumed  to  be  isotropic  and  can  be  measured  experimentally,  as  suggested  by 
existing  methods  [16].  For  the  analysis  of  the  unit  cell,  we  need  to  understand  the  relationship 
between  the  overall  loading/boundary  conditions  (OBC)  and  the  ligament  boundary  conditions 
(LBC).  The  LBC  varies  depending  on  the  location,  size  and  orientation  of  the  unit  cell  under  a 
certain  OBC. 

The  LBC  can  be  determined  by  the  following  bulk  analysis.  The  first  step  is  to 
generate  an  imaginary  cube  inside  which  a  unit  cell  of  the  foam  is  located,  as  depicted  in  Figure 
18.  The  imaginary  cube  can  be  generated  by  merging  a  cube  and  a  unit  cell  geometrically  by  the 
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Figure  18.  An  Imaginary  Cube  Inside  which  a  Unit  Cell  of  the  Foam  is  Located. 

ANSYS.  Overall  loading  and  boundary  conditions  are  then  applied  to  the  faces  of  the  bulk  cube. 
For  a  uniform  loading  condition,  one  of  the  faces  is  fixed  and  the  opposite  side  is  loaded, 
whereas  the  other  four  faces  are  free  to  deform.  Under  the  OBC,  the  bulk  cube  deforms 
isotropically  according  to  the  bulk  moduli  and  the  Poisson’s  ratio,  and  the  FEA  calculates  the 
deformations  and  stresses  point  by  point.  The  LBC  are  then  determined  by  the  results  on  the 
points  that  coincide  with  the  unit  cell  of  the  foam.  The  second  FEA  is  run  by  assigning  the 
recorded  displacements  to  the  unit  cell  as  the  LBC. 

The  current  method  is  used  for  the  moduli  back-calculation  of  the  carbon  foam. 
Effective  elastic  moduli  are  calculated  by  solving  six  different  cases  under  uniform  axial  and 
shear  strain  loadings  ( et).  For  each  uniform- strain  boundary  condition,  effective  stresses  (d,) 
are  calculated  by  taking  a  volumetric  average,  as  follows: 

J  a,.(x,  y,  z)  dxdydz 


a 


V 


,  (/  =  /,■••  ,6) 


(48) 


Components  of  6x6  effective  stiffness  matrix  ( [c (/  J )  and  effective  compliance  matrix  ( [sT  J)  are 
then  obtained  by  the  following  calculation: 
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where 


fej=[c, ,]{?,}  and  [sj=[cj' 
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Finally,  the  three-dimensional  effective  elastic  moduli  are  calculated  by  the 

following: 
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The  overall  effective  moduli  of  the  foam  are  then  calculated  by  repeating  the  above  procedure 
with  various  sizes  and  orientations  of  the  unit  cells  in  the  three-dimensional  space.  A  statistical 
approach  can  be  used  to  handle  a  huge  selection  of  random  sizes  and  orientations. 


2.3  SUMMARY 

The  emerging  ultra-lightweight  material,  carbon  foam,  is  modeled  with  the  three- 
dimensional  microstructures  to  develop  a  basic  understanding  of  the  performance  of  open-cell 
foam  materials.  The  model  will  describe  the  deformation  behavior  accurately  and  will  be  used  to 
investigate  the  failure  mechanism  of  the  cell  ligaments. 
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Because  of  the  randomness  and  complexity  of  the  microstructure  of  the  carbon  foam,  the 
representative  cell  ligaments  are  first  characterized  in  detail  at  the  microstructural  level.  The 
microstructural  characterization  will  then  be  correlated  with  the  macroscopic  bulk  properties  by  a 
statistical  approach.  A  series  of  databases  will  be  collected  for  various  sizes  and  spatial 
orientations  of  the  cell  ligaments,  as  well  as  the  property  variation  due  to  the  graphic  alignment 
along  the  longitudinal  direction  of  the  ligaments. 

A  computer  program,  “3D  Foam,”  will  be  developed  to  predict  longitudinal  and 
transverse  displacements  as  well  as  rotations,  and  calculate  the  reliable  stress  and  strain 
distributions  along  the  ligaments.  Because  of  slenderness,  each  ligament  can  be  considered  as  a 
beam,  and  the  tetrahedral  cell  microstructure  with  four  ligaments  as  a  frame  structure.  The  four 
beams  are  located  in  three-dimensional  space  under  arbitrary  loading  conditions.  The  cross 
section  of  the  beam  varies  in  size  and  material  properties  in  the  longitudinal  direction  along  the 
ligament.  A  tool  integrating  the  process  model  along  with  the  micro-  and  macro-analysis  of  the 
carbon  foam  will  lead  to  an  optimal  process  design  to  improve  foam  quality  and  to  reduce  cost. 

2.4  CONCLUSIONS  AND  RECOMMENDATIONS 

An  analytic  model  was  developed  for  the  unit-cell  of  carbon  foam  utilizing  variable 
material  properties  and  ligament  cross-section  geometry,  which  are  consistent  with  an  open-cell 
foam.  The  model  will  be  used  to  correlate  the  microstructural  properties  such  as  graphitic 
alignment,  porosity  and  ligament  structure  with  bulk  properties  that  are  being  measured  by 
mechanical  tests.  Experimental  validation  of  the  model  has  yet  to  be  completed.  The  validation 
effort  will  include  measurement  of  bulk  materials  properties  and  observation  of  ligament 
deformation  using  the  miniature  test  fixture  in  the  SEM. 
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3.  MODEL  DEVELOPMENT  OF  CARBON  FOAM  BLOWING  PROCESS 


3.1  LITERATURE  REVIEW 

Considerable  literature  research  in  theoretical  and  experimental  areas  was  conducted  on 
polymer  foam  formation  processes  [17-22].  Polymer  foam  can  be  produced  by  a  wide  variety  of 
processes:  expandable  beads,  injection  molding,  and  extrusion.  However,  all  of  these  processes 
have  one  basic  phenomenon  in  common:  the  nucleation  and  subsequent  nonisothermal  growth 
of  bubbles  upon  sudden  supers aturation  of  a  solution  consisting  of  a  gas  dissolved  in  the  melt 
polymer.  This  is  similar  to  the  carbon  foam  forming  process.  In  the  carbon  foam  blowing 
process,  a  gas  dissolves  or  is  trapped  in  the  pitch  under  high  pressure.  As  the  high  pressure 
releases,  the  bubbles  in  the  supersaturated  solution  form,  grow,  and  coalesce.  The  bubble  walls 
between  cells  open  up,  which  results  in  forming  the  open  cells  [23].  The  current  models  of 
polymer  foam  processes  are  based  on  the  following  assumptions: 

(1)  Nucleation 

The  nucleation  in  polymer  foaming  processes  can  be  classified  as  two  types: 
homogeneous  and  heterogeneous.  The  classical  nucleation  is  used  to  describe  the 
phenomena. 

(2)  Foam  bubble  growth 

(a)  The  bubble  is  spherical. 

(b)  The  gas  inside  the  bubble  follows  the  ideal  gas  law. 

(c)  The  gas  concentration  varies  only  with  the  radial  position  of  the  sphere  and  time. 

The  relationship  between  the  gas  pressure  in  the  bubble  and  the  residual  gas 
concentration  on  the  liquid  layer  surrounding  the  bubble  follows  Henry’s  law. 

(d)  The  melt  polymer  properties  are  independent  of  the  gas  concentration. 
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(e)  Because  the  process  is  so  quick,  it  can  be  considered  as  an  isothermal  process. 

(f)  The  inertial  effects  are  negligible. 

(g)  The  melt  polymer  is  considered  as  a  Newtonian  fluid. 

(3)  Mold  filling 

In  the  mold  filling  process,  the  mass  and  momentum  conservation  equations  are  used  to 
simulate  the  process.  As  the  thickness  of  a  foam  is  usually  much  smaller  than  the  length 
and  width,  the  Hele-Shaw  equation  is  often  used  to  simplify  the  three-dimensional  flow. 
As  the  bubbles  grow,  the  foam  volume  increases  (density  of  foam  reducing)  to  fill  the 
mold  cavity. 

3.2  RESEARCH  APPROACH 

Initially,  homogeneous  nucleation  will  be  assumed  in  the  model.  The  classical  nucleation 
theory  will  be  used  in  the  model  [24,25]. 

Most  work  in  this  research  will  focus  on  the  bubble  growth,  because  the  bubble  shape, 
size  and  distribution  determine  the  foam  properties  [25-29].  The  current  models  in  the  polymer 
foam  process  may  be  used  for  the  initial  bubble  growth  after  modifying.  However,  they  cannot 
be  used  for  the  carbon  foam  process  that  is  an  open-cell  foam  forming  process.  After  the  bubbles 
touch  each  other,  the  walls  of  the  impinging  bubble  collapse.  This  forms  the  open-cell  foam.  In 
the  fluid  mechanics  theory,  the  determination  of  the  cell  structure  is  a  free-surface  problem.  This 
is  quite  a  challenging  topic  because  the  location  of  the  free  surface  is  not  known  a  priori ,  and  the 
shape  of  the  free  surface  influences  the  flow  through  a  complicated  nonlinear  boundary 
condition,  the  normal-stress  condition.  When  motion  is  steady  or  quasi-steady  at  a  free  surface 
between  liquid  and  gas,  the  boundary  conditions  are  as  follows:  (a)  no  flow  normal  to  the 
surface,  (b)  no  tangential  stress  on  the  surface  if  the  surface  tension  gradients  are  negligible,  and 
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(c)  balancing  of  gas  pressure,  liquid  pressure  and  normal  viscous  stress  of  the  liquid  with  the 
capillary  pressure  that  is  the  product  of  the  surface  tension  and  mean  curvature  of  the  surface. 

All  iterative  schemes  employ  a  similar  strategy.  First  a  location  of  the  free  surface  is  chosen, 
either  by  an  informed  guess  or  on  the  basis  of  the  previous  iterations.  The  governing  equations 
of  mass  and  momentum  are  solved  for  the  velocity  and  pressure  fields  in  the  liquid,  but  only  two 
of  the  three  boundary  conditions  are  satisfied.  The  residual  in  the  third  boundary  condition  is 
used  to  decide  how  to  alter  the  location  of  the  free  surface.  The  calculation  process  is  repeated 
until  the  calculation  error  is  smaller  than  a  criterion.  This  scheme  is  complex  and  time- 
consuming  because  meshes  have  to  be  regenerated  in  the  calculated  domain  at  each  calculation 
step  (one  time  step  consists  of  many  calculation  steps),  and  convergence  is  often  difficult  to 
obtain.  On  the  other  hand,  the  model  used  in  polymer  foam  forming  is  much  simpler  than  the 
free-surface  calculation.  Although  it  may  not  be  as  accurate  as  the  previous  one,  it  can  provide 
basic  information  about  the  foam  growth  and  process  parameters.  At  this  stage  a  similar  method 
in  polymer  foam-forming  is  used  after  modification.  At  the  initial  bubble  growth,  the 
assumption  of  bubble  spherical  shape  is  acceptable,  as  the  following  calculation  indicates: 


dR  _  (Pg  ~pl)R  o 
dt  4p  2p 


d  (  4;t  PgR3  " 

dt  3  9tT 

V  7 


47tR2D— 

5r 


lr=R 


dc  dR  R2  dc  _  1  d  2  dc 
dt  dt  r2  dr  r2  dr  dr 


49 


The  initial  and  boundary  conditions  are  as  follows: 

R(0)=Ri 

Pg(0)=Pg,i 

c(r,0)=Ci(r) 
c(R,t)=KHPg(t) 
c(°°,  t)=c0 

where  c  is  the  concentration  of  the  dissolved  gas  in  melt  pitch,  D  is  the  diffusion  coefficient,  R  is 
the  radius  of  the  bubble,  r  is  the  radial  distance  from  the  bubble  center,  91  is  the  ideal  gas 
constant,  Pg  is  the  gas  pressure  inside  the  bubble,  Pl  is  the  pressure  of  the  melt  pitch,  and  o  is  the 
surface  tension  of  the  melt  pitch. 

After  the  bubbles  touch  and  collapse,  the  equations  above  have  to  be  modified.  As  the 
spherical  shape  assumption  is  still  used,  R  still  is  the  radius  of  the  assumed  sphere,  and  the 
bubble  volume  will  be  the  segment  of  the  assumed  sphere.  The  gas  diffusion  area  will  be  the 
segment  area. 

The  mold-filling  process  of  the  carbon  foam  may  be  considered  as  an  isothermal  process 
because  the  process  takes  a  very  short  time.  The  mass  conservation  equation  is  as  follows: 

%!L+V(P«llV)=0 

where  pceii  is  the  density  of  a  cell  that  consists  of  the  bubble  and  the  melt  pitch  surrounding  it. 

If  the  mold  thickness  is  smaller  than  its  length  and  width,  the  Hagen-Poiseuille  (H-P)  equation 
can  be  used.  Combining  the  mass  balance  and  H-P  equation,  a  Poisson  equation  can  be  obtained. 


50 


A  finite  element  code  is  being  developed.  Two  element  formats  were  selected: 
tetrahedral  and  isoparametric.  The  tetrahedral  element  is  used  to  connect  the  analysis  to  the 
mechanical  property  of  the  foam.  The  isoparametric  element  will  be  used  in  the  future.  The 
subroutines  to  calculate  the  coefficient  matrix  have  been  finished. 

3.3  CONCLUSIONS 

The  model  used  in  simulating  the  polymer  foam  process  is  modified  to  handle  the  open¬ 
cell  case  in  the  carbon  form  process.  Nucleation,  bubble  growth  and  mold-cavity-filling  will  be 
included  in  the  model.  At  this  stage  the  spheretic  bubble  assumption  is  used.  Several 
subroutines  to  calculate  the  coefficient  matrix  of  FEA  have  been  completed. 
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LIST  OF  ACRONYMS 


Acronym 

Description 

EAM 

element  array  model 

FEA 

finite  element  analysis 

FPF 

first-ply-failure 

H-P 

Hagen-Poiseuille 

FBC 

ligament  boundary  conditions 

FPF 

last-ply-failure 

MKM 

modified  Kabelka’s  model 

MMPM 

modified  mosaic  parallel  model 

OBC 

overall  loading/boundary  conditions 

RVE 

representative  volume  element 

SAM 

slice  array  model 

SEM 

scanning  electron  microscope 
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